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ABSTRACT 

We calculate the long-term evolution of angular momentum in double white dwarf binaries 
undergoing direct impact accretion over a broad range of parameter space. We allow the rotation 
rate of both components to vary, and account for the exchange of angular momentum between 
the spins of the white dwarfs and the orbit, while conserving the total angular momentum. We 
include gravitational, tidal, and mass transfer effects in the orbital evolution, and allow the Roche 
radius of the donor star to vary with both the stellar mass and the rotation rate. We examine 
the long-term stability of these systems, focusing in particular on those systems that may be 
progenitors of AM CVn or Type la Supernovae. We find that our analysis yields an increase 
in the predicted number of stable systems compared to that in previous studies. Additionally, 
we find that by properly accounting for the effects of asynchronism between the donor and the 
orbit on the Roche-lobe size, we eliminate oscillations in the orbital parameters which are found 
in previous studies. Removing these oscillations can reduce the peak mass transfer rate in some 
systems, keeping them from entering an unstable mass transfer phase. 

Subject headings: Celestial mechanics, Stars: Binaries: Close, Stars: Mass Loss, Accretion, Methods: 
Numerical 
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1. Introduction 

Binary systems containing compact objects are 
of great interest in a variety of areas in astro¬ 
physics. Of particular importance are double 
white dwarf (DWD) binary systems where both 
components are white dwarfs, which may make up 


the largest fraction of close binary stars (Marsh et 
al. ||1995[ ). 


Following common envelope evolution, DWD 
binaries may emerge with a sufficiently small 
semi-major axis, allowing gravitational radiation 
to drive the stars closer together on an astro- 
physically interesting time scale. These close 
DWD systems are of prime importance as low- 


frequency gravitational-wave sources (Hils, Ben- 


der and Webbink 1990; |Hils and Bender 2000 


Nelemans et al. 2001) as well as progenitors of 
Type la Supernovae ( Maoz et al. ||2014 ). 

As energy loss due to gravitational waves drives 
the degenerate components of a DWD binary to¬ 
gether, it is possible for the system to enter into a 
stable semi-detached state, in which the less mas¬ 
sive component will fill its Roche lobe and be¬ 
gin transferring matter to its companion. Such 
systems may result in the formation of AM CYn 


(Nather et al. 1981; 

Tutukov and Yungelson 1996 

Nelemans et al. 2001 

), which have extremely short 


orbital periods. During mass transfer, a stream 
of matter is pulled from the donor star through 
the inner Lagrangian point. If the matter stream 
does not impact the surface of the companion star, 
then the mass lost from the donor is expected to 


settle into a disc (Frank, King and Raine 2002). 


Torques exerted between the discs and the compo¬ 
nent white dwarfs allow angular momentum stored 
in the disc to be transferred back to the orbit 


(Soberman et al. 1997 Frank, King and Raine 


2002). This allows matter in the disc to fall onto 


the companion star while slowing the orbit con¬ 
traction, increasing the potential for a stable, long- 
lived binary system. Q 

For some DWD binaries, the mass transfer 
stream will directly impact the surface of the com- 


1 We note that we have not considered the possibility of 
the white-dwarf components having some residual, non¬ 
degenerate, hydrogen-rich outer layers. In a recent anal- 
ysis, [S hen 1 12015[| showed that the inclusion of such effects 
and associated nova-like outbursts may predominantly lead 
to mergers for both direct-impact and disc accretion. 


panion star. In doing so, there is no longer an 
obvious mechanism to return the orbital angular 
momentum from the transferred mass to the orbit. 
The division between the disc accretion and direct 


impact has been studied before (see Sepinsky et al. 
2014 Marsh et al. |2004 Gokhale et al. |2007 and 

references therein, hereafter Paper I, Marsh2004, 
and Gokhale2007, respectively). 

As in Marsh2004 and Gokhale2007, we showed 
in Paper I that mass transfer through direct im¬ 
pact can either increase or decrease the semi-major 
axis of the system, depending on the amount of an¬ 
gular momentum and mass exchanged between the 
components. If direct impact drives the system 
apart, both Marsh2004 and Gokhale2007 showed 
that a stabilizing accretion disc is likely to be 
created. If direct-impact mass transfer decreases 
the semi-major axis, the mass transfer rate may 
eventually become unstable, which can result in a 
merger. If the total mass of the system is in excess 
of the Chandrasekhar limit, such merger events 


could lead to a Type la supernova (Woosley and 
Weaver |[l 986 1. 


Marsh2004 and Goklrale2007 both calculated 
of the long term evolution of such systems. 
Marsh2004 concluded that the population of DWD 
binaries is likely lower than previously anticipated 
as a result of the high percentage of systems un¬ 
dergoing direct-impact accretion that become un¬ 
stable. Gokhale2007 improved upon the analysis 
of Marsh2004 by permitting the spin of the donor 
to vary and by including the effects of tidal forces 
from the donor star (Marsh2004 only examined 
tidal forces arising from asynchronicity between 
the accretor and the orbit. In their calculation, 
the spin of the donor was fixed to that of the 
orbit throughout so that no tides arose from asyn¬ 
chronicity between the donor and orbit.) The 
modifications of Gokhale2007 resulted in an in¬ 
crease in the number of stable systems which can 
be seen in Figures 2 and 3 in section 4. In Paper 
I we built upon both of these studies, providing 
direct ballistic integrations of the mass transfer 
stream (as opposed to the previous method us¬ 
ing an approximation adopted from |Verbunt and| 
Rappaport (1988)). We showed that removing 
this approximation which accounts for the com¬ 
plex three-body dynamics of the ejected mass can 
have a significant impact upon the angular mo¬ 
mentum exchange. 
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In this paper, we apply the results of Paper I to 
determine the long-term evolution of double white 
dwarf binaries in a fully self-consisent manner. In 
section 2, we introduce the equations which govern 
the long-term evolution of DWD systems and dis¬ 
cuss the differences between this method and the 
method utilized by previous studies. In section 3 
we discuss some details of our long-term numerical 
integrations. In section 4 we discuss the results of 
our solutions and analyze the results in compar¬ 
ison to those of previous works. We conclude in 
section 5. 


2. Equations of Long-term Evolution 


2.1. Basic Assumptions 

Following Paper I, we consider a close binary 
system of two white dwarfs with masses Ma and 
Md, volume-equivalent radii of Ra and Rd, and 
uniform rotation rates f 1a and Qd with axes per¬ 
pendicular to the orbital plane for the accretor 
and donor, respectively. [^] We assume the mass 
of each star is distributed spherically symmet¬ 
rically. As in the analyses of Marsh2004 and 
Gokhale2007, we assume the binary to remain in 
a circular Keplarian orbit throughout its evolu¬ 
tion. The radius of each object is assigned follow¬ 
ing Eggleton’s zero-temperature mass-radius re¬ 


lation (equation 15 of Verbunt and Rappaport 
(1988)). We assume that both the donor and ac¬ 
cretor initially rotate sychronously with the orbit 
(Marsh2004, Goklrale2007)[^] We choose the initial 
semi-major axis of the orbit such that the volume- 
equivalent radius of the donor (Rd) is equal to the 
volume-equivalent radius of its Roche lobe as fit by 


Eggleton (1983). 


2 Throughout this paper, the subscripts “A” and “D” will 
correspond to the accretor and donor, respectively. 

3 We note that there may exist a substantial population of 
DWD with non-zero eccentricity and asynchronous compo¬ 
nents prior to the onset of direct impact accretion | |WilIems| 
|et al. ( 2007| ); |Bours et al. ] ( |2014[ > and references therein). 
For this analysis, we limit our scope to systems that are 
initialy circular. Initially circular systems are expected to 
remain circular due to tidal circularization. A more thor¬ 
ough treatment of the initial parameter-space including ec¬ 
centricity and asynchronicity will be presented in a forth¬ 
coming paper. 


2.2. Evolution of angular momentum 

As described in Marsh2004 and Gokhale2007, 
the orbital evolution of a circular binary system 
can be described by its orbital angular momentum. 
We begin by presenting each of the components 
that can lead to a change in the orbital angular 
momentum. 

The total angular momentum, J to t , of a binary 
system is given by the sum of the orbital angular 
momentum, J or b, and the spin angular momenta, 
J spin,A and J sp in,D , of the accretor and donor, re¬ 
spectively: 

Jtot — Jorb T Jspin.A T J spin, I) • (1) 


The orbital angular momentum is given by 


Jorb — 



MaMd 


( 2 ) 


where M = Ma + Md, G is the gravitational con¬ 
stant, and a is the semi-major axis of the system. 
The spin angular momenta of the accretor and 
donor are J sp in,A = kAM A R 2 A ^A and J sp in,D = 
kuMDR 2 D^D, respectively, where kA and ku are 
dimensionless constants depending upon the in¬ 
ternal structure of the accretor and donor, respec¬ 
tively. It follows that: 


Jt, 



MaMd 


+ kAM A R 2 A^A + kDMDR^D^D- 


(3) 


There are three effects that change the or¬ 
bital angular momentum over time: mass transfer 
(MT), tides, and gravitational radiation (GR). As¬ 
suming each effect is indepedent of the others, we 
can then write the total change of the orbital an¬ 
gular momentum as the sum of the changes due to 
each of the above effects, with subscripts, as noted 
above, for each of the respective components: 

Jorb Jorb,MT T Jorb,tides + Jorb,GR- (4) 

To determine the total change in the angular mo¬ 
mentum, then, we simply need to write the change 
due to each of the above effects. 

The change in orbital angular momentum due 
to GR for a circular orbit is given by: 


• 32 G 3 M a M d M t 

Jorb,GR — ~ 5 4 Jorb 

DC CL 


(5) 
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(see for example, Landau and Lifshitz 1975). 

Prior to the onset of mass transfer, the semi¬ 
major axis of the binary will have been shrinking 
due to the effects of GR as seen in equation ([5|. 
During this time, tidal coupling will act to circu¬ 
larize the binary as well as synchronize the spins 
of the stars with the orbit. At the onset of mass 
transfer, we assume the spins and orbit to be syn¬ 
chronized. As mass transfer begins, angular mo¬ 
mentum will be exchanged between the spins of 
the component stars and the orbit (see Paper I). 
Any resulting asynchronization between the spins 
and the orbit leads to tidal coupling, the strength 
of which will determine how much of the spin an¬ 
gular momentum of the accretor is returned to the 
orbit. This will ultimately affect the stability of 
the mass transfer process. 

As in Gokhale2007, we model the change in the 
binary orbital angular momentum due to tides by: 


j _ UaMaR 2 a ! kpM D R 2 D 

J orb,tides — ^ A i 

ta t d 


wd• (6) 


The first term on the right-hand side of this equa¬ 
tion represents the torque due to dissipative cou¬ 
pling upon the accretor, and the second term the 
torque upon the donor. These torques are parame¬ 
terized in terms of the synchronization time-scales 
of the accretor, ta, and donor, Tp and are linearly 
proportional to the difference between the compo¬ 
nent and orbit spin, 

omegai = Q,; — f i or b, with i £ {A, D} for the ac¬ 
cretor and donor, respectively. Here, fl or b is the 
angular velocity of the circular orbit: 


n 


orb — 



(7) 


and k,M,R'j are the moments of inertia of each 
component, with fcj being the inertial constant. 


The values of ta and Tp determine the strength 
of tidal coupling relative to MT and GR. As in 
Gokhale2007, we examine both the cases of ta = 
td = 10 15 years (weak tidal coupling) and ta = 
t d = 10 years (strong tidal coupling). The tidal- 
synchronization timescales are discussed further in 

o 


To determine J 0 rb,MT , we follow Paper I, which 
uses the ballistic mass transfer calculations of 


Sepinsky et al. (2010) to determine the instan¬ 


taneous effect of mass transfer on DWD systems. 


This method uses a fully self-consistent, conser¬ 
vative, ballistic model of the transferred mass to 
determine the orbital parameters of the system af¬ 
ter a single mass-transfer event. As seen in that 
paper, provided the mass transferred is small com¬ 
pared to the total mass of the binary system, the 
change in the orbital parameters per unit mass 
that results is directly proportional to the mass 
transfer rate: 


Jorb,MT 


Ad orb .b 
Mp 


Mp 


( 8 ) 


where A J or b,b is the change in the orbital an¬ 
gular momentum for the DWD as calculated by 
the above ballistic model for a single mass trans¬ 
fer event ejecting a particle of mass Mp. In this 
method, changes in J or b,MT (and changes in all or¬ 
bital parameters) are calculated at each time-step 
by integrating the three-body system consisting 
of the two stars and the discrete particle repre¬ 
senting the mass transfer stream. The change in 
the J 0 rb,MT per unit mass transferred is indepen¬ 
dent of the mass of the ejected particle as long 
as Mp « Md,Ma■ For the calculations here, 
we use Mp = 10 _8 M 0 . The rate of change of 
Jorb.MT is then determined by multiplying by the 
current mass-transfer rate of our evolving DWD, 
Mp- The calculation of Mp will be discussed in 
section 12.51 


2.3. Differential equations for long-term 
evolution 


Using the above rates of change for the orbital 
angular momentum, we are now ready to develop 
the equations for long-term evolution. It follows 
from equation ([2]) that 



/, x m d 1 a 
[1 - q) M- D + 2a 


(9) 


for conservative mass transfer (M = 0)|^] We let 
q = Mp/Ma- 


2.3.1. Evolution of the semi-major axis 

As in section [2~2| for the orbital angular momen¬ 
tum, we examine the changes of the semi-major 
axis due to mass transfer, tides, and GR. We as¬ 
sume that each effect is independent, and write the 

4 Recall that we assume the orbit remains circular through¬ 
out. 
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total change in the semi-major axis due to each of 
the above effects, respectively, as: 


a = cimt + audes + a-GR- 


( 10 ) 


To calculate omt, we utilize the model for mass 


transfer developed in Sepinsky et al. (2010). The 


differences between this model and models used in 
other analyses will be discussed in the section [274] 
Using this method, we calculate the time rate of 
change of a due to mass transfer as 


Aab 

aMT ~ Wp D ' 


( 11 ) 


Here, Actf, is the change in the semi-major axis 
for the DWD as calculated by the above ballistic 
model of Sepinsky et al. (2010) as described in 


section m for the orbital angular momentum. 

Next we calculate dudes and clgr ■ We cannot 
use the standard tidal prescription for changes in 


the semi-major axis following Hut (1981) because, 


following Marsh2004 and Gokhale2007, we used 
a different metric for changing the spin angular 
momentum due to tides (equation [6] ). Instead, we 
must develop dudes from that angular momentum 
change. Holding the donor mass constant ( Md = 
0) for the final two terms of equation (10) we can 
combine equation Q with equation (|9 ) to obtain: 


Jorb , 


tides 


J orb, 


GR 


J, 


orb 


J, 


orb 


n O'tides,GR 

2a 


( 12 ) 


where dudes,GR is the total change of the semi¬ 
major axis due to the combined effects of tides 
and GR. Assuming as before that changes to the 
semi-major axis due to tides and changes due to 
GR are independent we can re-write the above as: 

Jorb, tides Jorb,GR 1 . 1 . /iq\ 

- J -1- Y~ - _ O ~ atid es + 7Z~ a GR- U'jJ 

t/ or b a orb zci 

Using equations 0,0 , and ([6]) , it follows that: 


dtides — 



( ^aMaRa 

X - UJ A 

\ ta 


koMoRp 

td 


LOd 


(14) 


and 


&GR = 


64 G 3 M a M d M 
5 c 5 a 3 


(15) 


For our purposes, we rewrite the two tidal com¬ 
ponents in terms of the rotation rates of the accre- 
tor and the donor relative to the orbital angular 
velocity, which we define as f a and fp>, respec¬ 
tively. The rotation rates can be written in terms 
of the rotational angular velocities of the stars, fl; 
and the angular velocity of the circular orbit, fl or b' 


h - 1 = 


^orb i 


n 


orb 


n 


orb 


(16) 


Using equations 0 and ( fl6| ), we can rewrite 
equation (14) as 


2 M 


&tides — 


where 


and 


ola = 


old = 


clMaMd 

k aMaR"a 

ta 

k d M d R- 


d 


td 


(«A + Old) 

(17) 

if A- 1) 

(18) 

ifn-1). 

(19) 


Finally, we can insert equations (11), (15), and 
into equation (JToJ) to obtain the equation for 
the evolution of the semi-major axis with time: 


a = 


A at 
Mp 


M d + 


2 M 

clMaMd 


{OLA + 0l D ) 


64 G 3 M a M d M 
5 c 5 a 3 


( 20 ) 


2.3.2. Evolution of the component rotation rates 


Next, we find the equations for the evolution of 
the component spins, f a and fa- Like the changes 
to the semi-major axis (equation [lO] ), the changes 
in Ja and /d can be separated into three compo¬ 
nents: 

fi — fi,MT T fi,tides 4“ fi,GR (^1) 


with i £ { A , D}. 

Analogous to omt in equation the change 
in fi due to mass transfer, fi,MT can be written 
as: 


fi, 


MT — 


Afbd 

Mp 




D- 


( 22 ) 


Here A fb,i is the change in the rotation rates of 
star i resulting from a single mass transfer event 
in the formulation of Sepinsky et al. (2010) as de¬ 
scribed in section El 
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As in equation ([3]), the spin angular momentum 
of each component can be written as 


Jspin,i — U Adj R, f f / 1 1 


(23) 


orb • 


where we have substituted f t from equation (|16|. 

Since we have already determined the change in 
fi due to mass transfer (equation [22] ), we can de¬ 
termine fi,tides and fi,GR by differentiating equa¬ 
tion (23) with the mass held constant: 




spin,', 


Mi 


— kiAfiR^ Qorb(fi,tides fi,GR) 


3 1 

2 kiMiRi fl or bfi — {(itides o>gr) ( 24 ) 

where the second term arises due to the de¬ 
pendence of fl or b on the semi-major axis (equa¬ 
tion [7]). We note that, because we are holding the 
mass constant, the second term depends only upon 
changes due to tides and GR, and not changes due 
to mass transfer. Changes to f, due to the effect 
of mass transfer on the semi-major axis are fully 


accounted for by equation (22). 


Since we do not include any GR effects on the 
spin angular momentum of the components, con¬ 
servation of angular momentum dictates that any 
changes in the spin angular momentum of a com¬ 
ponent must be equal and opposite to the changes 
in the orbital angular momentum of the system 
due to tides acting on that component. Combin¬ 
ing equations ©. 0- and (|16[), we have: 


J. 


spin,i 


k t M 2 Rf 


M, 


Vorbifi ~ !)• (25) 


By combining equations (24) and (25) and re¬ 
arranging, we can write fcR,i + fades,i as: 


fGR,i "F ftides,i — 


(fi - 1) 


3 h 

2 a 


2 Ad 

~P + “ 71 r nr ( aA + a D) 
a MaMd 


(26) 


tions for the evolution of /a and fp>\ 


U = A h A M D - fA 1 


Adp 


ta 


3 fA 
2 a 


-P+- 


M 


cl 


(a a+ an) 


3/d 

2 a 


fo = 

-P + 


Afb,P 

Alp 


Mp — 


/d-1 

td 


2 M 
a AIaMd 


(a a + cup) 


(28) 


(29) 


2.3.3. Evolution of the eccentricity 

We assume that the eccentricity of these sys¬ 
tems is zero prior to the onset of direct-impact 
accretion. During the three-body integration of a 
single mass-transfer event, it is possible for the bi¬ 
nary to develop a small eccentricity. However, for 
simplicity and in accordance with Marsh2004 and 
Gokhale2007, we force the eccentricity to remain 
zero throughout. Because the system begins in a 
circular orbit, and due to the action of tidal forces 
and gravitational wave emission which both act 
to circularize the orbit, it is unlikely for any sig¬ 
nificant eccentricity to develop. In order to keep 
this orbital angular momentum in the orbit, we 
manually set our new semi-major axis to: 

a = a/( 1-e 2 ) (30) 


where e is the eccentricity gained by the system 
during mass transfer and a/ is the semi-major axis 
of the eccentric orbit. Using this modification, we 
force: 

e = 0, (31) 

and still conserve orbital angular momentum. A 
more thorough analysis in which we allow eccen¬ 
tricity to vary throughout the entire calculation 
will be presented in a forthcoming paper. 


2.3.4- Adass transfer rates 


where we have used equations (15), (17)- (19), and 
let: 


P = 


64 G 3 Ad A Ad D Ad 
5 c 5 a 3 


(27) 


Following the form of equation (21), we can 
combine equations (22) and (26) to write the equa- 


Finally, the changes in the masses of the accre- 
tor and donor are given, respectively, by: 


dAdA 

dt 

dAlp, 

dt 


-Md (32) 

= Md (33) 
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The mass loss rate of the donor, Md is obtained 
as described in section [2751 


Together, the set of equations (201, (281, (291, 


(31), (321, (33) can be integrated in time to calcu¬ 
late the evolution of the system. 


2.4. Differences from previous analyses 


There are several differences between our treat¬ 
ment of mass transfer and those of Marsh2004 
and Gokhale2007. In order to calculate the an¬ 
gular momentum exchange during mass transfer, 
the studies of Marsh2004 and Gokhale2007 both 
utilize a numerical prescription based on |Verbunt| 


and Rappaport (1988). In that formulation, it is 


assumed that the angular momentum transferred 
from the orbit to the spin of the accretor is exactly 
equal to the angular momentum of the ballistic 
particle in a circular orbit around the donor at 
its average radius during its motion from donor to 
accretor. Where the analysis of Marsh2004 keeps 
the spin of the donor fixed, Gokale2007 does allow 
the spin of the donor to vary, and notes that in 
doing so, the number of stable systems increases. 
However, the allowance of variation in donor spin 
is not done self-consistently. By introducing the 
spin angular momentum of the donor, there are 
now three separate sources/sinks of angular mo¬ 
mentum: the spin of the donor, the spin of the 
accretor, and the orbit. Paper I showed that an¬ 
gular momentum is transferred between each dur¬ 
ing mass transfer, and that the fraction of angular 
momentum transferred between each is strongly 
dependent on the ballistic trajectory of the trans¬ 
ferred mass, and hence the system properties. In 
this paper, we apply the ballistic calculations of 


Paper I to lift the dependence on the Verbunt and 


Rappaport (1988) approximations to more accu¬ 


rately determine the flow of angular momentum 
between the component spins and the orbit. 


Since orbital angular momentum changes are 
directly linked to changes in the mass ratio and 
semi-major axis, much can be learned from anal¬ 
yses such as that of Gokhale2007 and Marsh2004. 
However, a more complete analysis demands a 
thorough consideration of the spin angular mo¬ 
menta of the stars as well due to their indirect 
effect upon the system properties through tidal 
coupling. This is accomplished via the ballistic 
calculations described above, which evaluate the 


three-body problem throughout the evolution of 
the system to determine the precise effect mass 
transfer has on the evolution of the system. These 
calculations take into account not only the imme¬ 
diate feedback on the orbit and spin of the accretor 
during ejection, but also (1) the gravitational ef¬ 
fect on the orbit during the mass transfer process, 
(2) a calculation of the precise moment when the 
particle impacts the surface of the accretor, (3) 
the instantaneous properties of both the accretor 
and particle at impact, including the angle of im¬ 
pact, and (4) correctly divides the momentum of 
the impacting particle between linear and angular 
momentum based upon the angle of impact and 
the rotation rate of the accretor. We can then ex¬ 
amine the angular momentum exchange between 
all components, the spins of both stars and the 
orbit, and do so in a way that allows the spins of 
both stars to vary self-consistently. 

2.5. Calculation of mass-transfer rate 

As in Marsh2004, we define the overfill of the 
Roche lobe as: 


A = Rd — Rl 


(34) 


where Rd is the radius of the donor and Rl 
is the radius of the donor’s Roche lobe. The 
way in which the mass-transfer rate varies with 
A has been investigated in many analyses (see, 

(1972); 


for example, Paczynski and Sienkiewicz 


Webbink (1977); Savonije (1978)). In accordance 


with Marsh2004, we approximate the mass trans¬ 
fer as adiabatic (Webbink 1984| ). In the adiabatic 
regime the mass transfer rate is given by: 


87r 3 /5 Gm e \ 3 ^ 2 5 / 2 

M D = ( , q ) {ii e m n ) 5/2 x 


A J 


9 V h 2 
1 / 3/iMc>\ 3/2 


Porb \5 tlR 2/ yj 0 , 2^2 — 1) 

for A > 0 and zero for A < 

(1984); [Chandrasekhar 


(35) 


from 

Webbink 

Webbink 

(1977) 


0 (using results 
(1967); 


Here, 


is the mass of an 
electron, m n is the mass of a nucleon, /j, e is the 
mean number of nucleons per free electron in the 
outer layers of the donor (assumed here to be two), 
P or b is the orbital period, ry. = i?y/a, and /r and 
02 are given by the following: 

Md 


H = 


Mj\ + Md 


(36) 
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a 2 = 


1 - M 


L 'L1 


(l-*m ) 3 


(37) 


where Xli is the distance from the center of the 
donor to the inner Lagrangian point of the donor, 
in units of the semi-major axis (Webbink|l977). 


2.6. Calculation of Roche lobe 

In the case of Marsh2004, where the rotational 
velocity of the donor is fixed to the orbital velocity 
and where the orbit is circular throughout, the 
shape and volume of the Roche lobe depends only 
upon the mass ratio of the system. In that case, 
the approximation from |Egg1eton ; (1983), 


o 0-49g 2/3 , , 

L ’ Egs O. 65 2 / 3 + ln(l + g 1 / 3 ) ’ 

is sufficient. However, for eccentric and/or non- 
synchronous binaries, the shape and volume of the 
roche lobe also depend upon the eccentricity, the 
true anomaly, v, and the rotation rate, which can 
be expressed by the parameter: 


A(e,f,v) = 


fHl + e) 4 
(1 


(39) 


see 


Sepinsky et al. 


e cos v) 3 
(2007). Other analyses have 


considered the dependence of the Roche lobe on 
asynchronicity. In the case of a circular orbit, the 


Roche lobes described in Sepinsky et al. (2007) 


using the Roche-lobe parameter Ai directly re¬ 
flect the equipotential surfaces described in |Plavec~| 
(1958) and Kruszewski (1963). We note that this 


asynchronicity parameter Ai depends on the rota¬ 
tion rate of the component star in question, which 
means the size of the Roche lobe can be different 
for each object in the system. 

For our calculations, we perform a Monte Carlo 
integration to calculate the volume-equivalent 


Roche lobe radius (as outlined in Sepinsky et 
al. (2007)) over a two-dimensional grid in A-q 
parameter-space. We then use a bilinear inter¬ 
polation function to create a continuous function 
for the Roche lobe from the two-dimensional grid. 
The interpolated function is accurate to 0.01 per¬ 
cent or better in comparison to the Monte Carlo 
integration at equivalent points. 

In section |4.4| we examine the dependence of 
the results on the size of the Roche lobe. There, we 
compare the evolution of the DWD systems using 


the Eggleton function for the roche lobe, Rl, E gg , 
to the evolution including the asynchronicity pa¬ 
rameter which we name RL,asynch- 

2.7. Synchronization time-scale 

As the spin angular momentum of the accre- 
tor increases due to impact from the mass trans¬ 
fer stream, tidal coupling will redistribute some 
of this angular momentum back into the orbit. 
The strength of this tidal coupling is dependent 
upon the synchronization time-scales, ta and T 75 , 
as seen in equation (| 6 |. [Campbell | ([l984]) calcu¬ 
lated this timescale as: 


tc = 1.3 x 10 


7 (Ma \ 2 (_a_ 


V Mo) \Ra 


Ma/Mq 

La/Lq 


6 n .f rtr -| 5/7 

yr 
(40) 

Following Marsh2004, we retain the scaling 
with mass ratio and orbital separation, but, be¬ 
cause of the uncertainty of the overall magnitude 
of the time-scale, we allow the magnitude to vary. 
We define the overall magnitude by the time-scale 
at the moment of first contact and assume, as in 
Marsh2004 and Gokhale2007, that 

«<*>-<* (£)’(£)’ (41) 


and 


= < 42 > 


where Ca and Cjj are constants defined such that 
774 ( 0 ) = ta.i and T£>(0) = tdj ■ Here, taj and 
Tnj are the initial synchronization time-scales. 
As in Gokhale2007, we explore two different val¬ 
ues for the synchronization time-scale at contact: 
10 15 years (very weak tidal coupling) and 10 years 
(very strong tidal coupling). Much work on tidal 
synchronization in DWDs has been done (see, for 


, Valsecchi et al. (2013); 

Fuller and Lai 

Burkart, Quataert, and 

Arras (2014)). 


These analyses have shown that shorter synchro¬ 
nization time-scales may be a better approxima¬ 
tion for systems with short orbital periods, as we 
consider in this analysis, which may lend credence 
to our 10 year synchronization time-scale. For sim¬ 
plicity, we use the same initial time-scale for the 
donor and accretor (taj = tjjj) and refer to the 
initial synchronization time-scale for both stars as 
simply r. 

















































We have now determined the set of differen¬ 
tial equations which governs the evolution of the 
DWDs undergoing direct-impact accretion, in¬ 
cluding the methods for calculating the Roche 
lobe size, the mass transfer rate, and the synchro¬ 
nization timescales which determine the strength 
of tidal coupling. We now have all the steps in 
place to solve for the long-term evolution of the 
systems. 


3. Numerical Solutions 


We integrate equations (20), (28), (29), (31), 


(32), (33) using an 8tlr order Runge Kutta or¬ 


dinary differential equation solver (Galassi et al. 


2006). Excluding the losses due to gravitational 


radiation, the total energy and momentum of the 
system are conserved throughout the integration 
over the entire parameter space to one part in 
10 -3 or better. At the end of the evolution, we 
Calculate Jconserve ~ Jtotal Jtotal.I T Jtotal,GR 
where Jtotal.I is the initial total angular momen¬ 
tum and Jtotal ,gr is the total angular momentum 
lost to gravitational radiation. If angular momen¬ 


tum is perfectly conserved, we expect J, 
be equal to zero. 


conserve 


to 


Figure 1 shows the final value of Jconserve for 
both synchronization time-scales at contact over 
the entire parameter space of interest for solutions 
calculated using both the Eggleton approximation 
(left), Rl,E gg, and R L ,as V nch (right) for the calcu¬ 
lation of the Roche lobe. The colors correspond 
to different values of Jconserve as described in the 
caption. In general, systems with a smaller to¬ 
tal mass conserve total angular momentum better 
than systems with a higher mass. For many of 
the systems in the lower right, this is due to the 
fact that they reach a stable configuration (disc 
accretion) early in the integration. Therefore, the 
integration needs to run for only a few orbits. Sys¬ 
tems where Jconserve is larger tend to be systems 
where the numerical integration runs for many or¬ 
bits, allowing systematic errors in the numerical 
integration to accumulate. Even so, all systems 
presented in this paper conserve total angular mo¬ 
mentum to better than 1% throughout the entire 
evolution. 


We integrate over a period of 1 Gyr. As in 
Marsh2004, if Mjj exceeds 0.01 M 0 /yr at any 
point during the integration, the integration stops 
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Fig. 1.— Angular momentum conservation over entire 
parameter space for the two synchronization time-scales. 
The left-hand panels show the solutions obtained using 
Rl ,Egg for the calculation of the Roche lobe and the right- 
hand panel shows the solutions using RL.asynch • Calcu¬ 
lated here is Jconserve — Jtotal Jtotal, I 4 Jtotal, GR ''t the 
end of the integration. Blue systems conserve angular mo¬ 
mentum to an order of 10 —3 ; green systems to an order of 
10 -4 ; cyan systems to an order of 10 5 ; magenta systems 
to an order of 10 —6 ; yellow systems to an order of 10 7 ; 
and black systems to an order of 10 —8 or better. 



and the system is discarded as unstable. 

As a system evolves, it is possible to pass back 
and forth through phases of mass transfer and 
phases of no mass transfer as the semi-major axis 
and the two masses change. If the semi-major axis 
increases enough for mass transfer to stop alto¬ 
gether, the integration proceeds (with Mp = 0) 
until the action of gravitational radiation shrinks 
the orbit sufficiently for mass transfer to resume. 


As in Paper I, in this work we are only inter¬ 
ested in direct-impact mass transfer where the par¬ 
ticle impacts the surface of the accretor within one 
orbital period. In this case, the evolution of the 
orbital parameters is determined by the differen¬ 
tial equations presented in section 2. If the parti¬ 
cle does not accrete within one orbital period, it 
is likely that the accretion stream will eventually 
intersect with itself, ultimately leading to the for¬ 
mation of an accretion disc ( Sepinsky et al. |2010 1. 

Compared to direct-impact accretion, disc ac¬ 
cretion is known to provide a much more efficient 
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mechanism for redistributing spin angular momen¬ 
tum in the accretor back into the orbit through 
tidal coupling (Frank, King and Rai ne|[20 021. As 
a result, it is likely that once a system enters a 
phase of disc accretion it will remain in this phase 
or perhaps even become detached as the orbital 
separation continues to grow and the mass ratio 
decreases due to continued mass transfer. If at 
any point during our numerical integrations a disc 
is formed under these circumstances, the integra¬ 
tion stops and the system is assumed to be stable 
throughout its lifetime. A more thorough analy¬ 
sis which continues to track the evolution through 
potential disc phases will be performed in a future 
paper. This issue is discussed further in section 4. 


3.1. Maximum spin-rate of the accretor 

For weak tidal coupling, it is possible for the 
accretor to be spun up to its breakup rate (f Ik = 
■\JGMa/R\ )• We track f l a throughout our calcu¬ 
lations and note that only one system ever reaches 
this maximum spin for the accretor. This system, 
which has an initial donor mass of 0.1 M 0 and an 
initial accretor mass of 0.275 M 0 , is marked by a 
yellow dot in Figure 2. The accretor spin-rates 
in all other systems stay below the breakup rate 
throughout their evolution. 


3.2. Super-Eddington accretion 


As noted by Marsh2004, there are likely ranges 
of parameter space where systems, despite being 
stable, experience super-Eddington accretion at 
some point during their evolution. It is expected 
that a possible ultimate consequence of sustained 


super-Eddington accretion is a merger (Han and 
Webbink 1999 ?; Marsh et al. 2004), the same 


result as a dynamically unstable system, as they 
eventually can reach very high mass-transfer rates. 

As in Marsh2004, we calculate the Eddington 
accretion rate using a modified form of the calcu¬ 


lation used by Han and Webbink (1999): 



o . 

0 0.5 1 
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Fig. 2.— Results of evolution with initial tidal synchro¬ 
nization time-scale of r = 10 15 years using the Eggleton 
approximation for the Roche lobe. Red systems are stable 
throughout their entire evolution and have a total system 
mass which is sub-Chandrasekhar ( < 1.44 Mq). Green 
systems have a total mass which is sub-Chandrasekhar 
but went through a phase of super-Eddington accretion at 
some point during the evolution. Blue systems are stable 
throughout and have a total system mass which is super- 
Chandrasekhar. Magenta systems have a total mass that 
is super-Chandrasekhar and undergo a period of super- 
Eddington accretion. Black systems are unstable, mean¬ 
ing the mass-transfer rate exceeded 0.01 Mq/ yr. The solid 
yellow line illustrates the widely-used q > 2/3 instability 
boundary (see, for example, Marsh2004). The solid black 
line shows the boundary between disc and direct-impact 
accretion for initially synchronous and circular binaries. 
The dotted black line and dashed black line illustrate the 
boundary between stable systems and super-Eddington sys¬ 
tems for the t = 10 15 year timescale in Marsh2004 and 
Gokhale2007, respectively (see Figure 5 in Marsh2004 and 
Figure 7 in Gokhale2007). 


point of impact, both measured in the co-rotating 
frame of reference. 

In the Marsh2004 analysis, a set of impact ve¬ 
locities and locations were pre-computed and were 
interpolated during the calculations to calculate Vi 
and v w . Instead we are able to calculate Vi and 
v w explicitly as part of the three-body integration 
of the mass transfer event. We track the mass- 
transfer and flag systems that eventually exceed 
the maximum mass-transfer as unstable. 


8tt Gtyi v cMa 

Edd = 7“ " i 7 77 77T 

°TVPL1 ~ q>a~ 2 V i + 2( V i “ v w) ) 

(43) 

where ot is the Thomson cross-section of the elec¬ 
tron, uip is the mass of a proton, Vj is the im¬ 
pact velocity of the accreted particle, and v u is 
the spin-velocity of the accretor’s surface at the 


4. Results 

Using the techniques described above, we com¬ 
puted evolution over a grid in Ma,Md parameter 
space to determine the long-term stability of var¬ 
ious systems. The grid was computed for two dif¬ 
ferent tidal synchronization time-scales at contact: 
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10 15 years and 10 years. 

4.1. Evolution of systems using Eggleton 
Roche lobe 

We first present the results of our numerical 
integrations using the standard Eggleton approxi¬ 
mation for the calculation of the Roche lobe (equa¬ 
tion |38| ), which does not take into account asyn¬ 
chronicity between the donor and orbit. The 
Eggleton approximation was was used in the anal¬ 
yses of both Marsh2004 and Gokale2007. 

Figure 2 shows the end result of systems with 
an initial synchronization time-scale of t = 10 15 
years using Rl,E gg- As mentioned in section 3, 
we assume that once a disc is reached, the system 
will remain stable for the remainder of its evo¬ 
lution, therefore we stop the integration once a 
disc is formed. With the exception of the unsta¬ 
ble black systems (whose evolution stops when the 
mass transfer rate exceeds the maximum limit of 
0.01 Mq/yt), the evolution of all systems in this 
plot is stopped when the system reaches a phase 
of disc accretion. The time it takes to reach this 
phase of disc accretion varies from system to sys¬ 
tem. 

The solid black line in Figure 2 shows the 
boundary between disc and direct-impact accre¬ 
tion for initially synchronous and circular bina¬ 
ries. All systems to the right of and below this 
line begin evolution in a disc phase. Therefore, 
these systems are assumed stable from the onset 
of evolution and do not evolve past t = 0. Since 
these systems do not evolve in time, the angu¬ 
lar momentum is perfectly conserved, as shown by 

0.5 
0.4 

2 0,3 
sf 0.2 
0.1 

0 

0 0.5 1 

M a (M e ) 

Fig. 3.— As Figure 2, but with r = 10 years. 


the black systems in the lower right of the plots in 
Figure 1. All systems to the left of and above this 
line begin evolution in a direct-impact phase. We 
will discuss how we handle disc accretion further 
in section 4.2. 

Figure 3, analogous to Figure 2, shows the 
end results of systems with initial synchronization 
time-scales of r = 10 years using Rl,E gg for the 
Roche lobe calculation. 

The sub-Chandrasekhar, sub-Eddington sys¬ 
tems shown in red in Figures 2 and 3 are ex¬ 
pected to remain stable throughout their lifetimes 
and are therefore characterized as likely AM CVn 
progenitors. As observed from the dotted and 
dashed black lines in Figures 2 and 3, the region 
of parameter space occupied by such systems in 
Gokhale2007 (dashed line) increased the number 
of stable systems compared to Marsli2004 (dot¬ 
ted line); our analysis reveals a further increase 
in the extent of this region of parameter space. 
This observation holds true for both the 10 15 year 
time-scale (Figure 2) and the 10 year time-scale 
(Figure 3). 

Gokhale2007 noted that by allowing the spin of 
the donor to vary and by including the resulting 
effects of tidal coupling between the donor’s spin 
and the orbit, the number of stable systems in¬ 
creased compared to the analysis of Marsh2004 (as 
observed by comparing the dotted versus dashed 
lines in both Figures 2 and 3). As described in sec¬ 
tion 2.4, the main difference between our analysis 
and that of Gokhale2007 is in the treatment of the 
way we handle angular momentum evolution dur¬ 
ing the direct impact accretion process. In observ¬ 
ing the increase in the number of stable systems 
compared to the analysis of Gokhale2007, we con¬ 
clude that by utilizing a mass transfer treatment 
that allows the rotation rates of both components 
to vary and self-consistently accounts for the ex¬ 
change of angular momentum between the spins 
of the components and the orbit, we are able to 
increase the number of stable systems. 

Figures 2 and 3, produced using Rl,E gg for the 
Roche lobe calculation, isolate the effect of our 
different treatment of the mass transfer process 
has upon the stability of systems in comparison 
with Marsh2004 and Gokhale2007. We will discuss 
the differences resulting from using RL,asynch in 
section 4.4. 
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Forming a boundary between the stable (red 
and blue) and unstable (black) systems are the 
green and magenta systems, which experience 
super-Eddington accretion at some point dur¬ 
ing the integration. Unlike the black systems, 
whose mass transfer exceeds the allowed limit 
of O.OlM 0 /yr and are therefore categorized as 
unstable, these green/magenta super-Eddington 
systems are categorized as stable. However, it is 
likely that systemic mass loss is needed for a bi¬ 
nary to survive phases of super-Eddington accre¬ 
tion. Han and Webbink (1999) argues that such 
mass loss may lead to a common envelope sur¬ 
rounding the binary, which will ultimately lead to 
a merger. In this case, systems which experience 
super-Eddington accretion would be considered 
unstable systems. The analyses of Marsh2004 and 
Gokhale2007 both note that super-Eddington sys¬ 
tems are expected to be unstable. We make the 
same assumption but note that a more thorough 
treatment of the physics governing systems as they 
pass through any phases of super-Eddington accre¬ 
tion is necessary to state with certainty whether or 
not such systems are ultimately stable or unstable. 


The green super-Eddington systems have a to¬ 
tal mass which is sub-Chandrasekhar, and the 
magenta super-Eddington systems are super- 
Chandrasekhar. If we assume super-Eddington 
systems are ultimately unstable, the magenta sys¬ 
tems can be categorized as possible Type la su¬ 
pernova progenitors. 

As the tidal synchronization time-scale is re¬ 
duced from r = 10 15 years in Figure 2 to r = 10 
years in Figure 3, we see the parameter space oc¬ 
cupied by stable systems grows (red). This is to 
be expected and in agreement with the analyses 
of Marsh2004 and Gokhale2007. Stronger tidal 
coupling will allow more spin angular momentum 
from the accretor to be transferred back into the 
orbit, which increases the semi-major axis, causing 
the mass-transfer rate to decrease. This results in 
an increase the stability of the systems in general. 


4.2. Disc accretion 

In both Figures 2 and 3 we see a portion of sta¬ 
ble, super-Chandrasekhar, sub-Eddington systems 
(blue) in the top right corners of our plots. In the 
analyses of both Marsh2004 and Gokhale2007, this 
portion of parameter space is occupied by super- 
Eddington systems, which would be magenta cir- 



0 50 100 


time (years) 


Fig. 4.— Various system properties for the first 150 
years of evolution of a system with = 0.38 Mq and 
Md = 0.8 Mq. Evolution using an initial synchronization 
time-scale of r = 10 15 years is shown in black; this system 
is unstable (see Figure 3). The same system with r = 10 
years is shown in green; this system is super-Eddington. 
t = 0.1 year for this system is shown in red; this system is 
stable. 


cles in our plots. Notice that the dotted and 
dashed lines of Marsh2004 and Gokhale2007 de- 
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crease as we move right across the plot, whereas 
our stability boundary is more extended. This is 
due to the fact that we do not follow the evolu¬ 
tion of the systems past the development of disc 
accretion, while Marsh2004 and Gokhale2007 both 
continue to track the evolution of systems through 
phases of disc accretion. 

Here we acknowledge the main limitation of our 
decision to assume that systems are stable once 
they reach a phase of disc accretion. The blue 
systems in the region of parameter space of inter¬ 
est here begin the integration in a disc phase, and 
are therefore labeled as stable. However, it is fea¬ 
sible that if these high donor-mass systems were 
allowed to evolve through the initial disc phase, 
they would eventually become super-Eddington 
and potentially unstable. In this case, this region 
of parameter space would more closely conform 
with the shape of Marsh2004’s and Gokhale2007’s 
analyses, both of which allow the systems to con¬ 
tinue evolving through any phases of disc accre¬ 
tion. As we noted earlier, a more thorough anal¬ 
ysis of this region of parameter space which in¬ 
cludes treatment of disc accretion will be studied 
in a forthcoming paper. 

4.3. Comparison of system evolution with 
different tidal time-scales 

Figure 4 shows the evolution in time of a sys¬ 
tem with initial masses of Ma = 0.8 Mg and 
Mjj = 0.38 Mq for the two different synchroniza¬ 
tion timescales at contact as well as a third much 
stronger tidal timescale, r = 0.1 year, for further 
illustration. These solutions are obtained using 
Rl, Egg to calculate the Roche lobe. For r = 10 15 
years, the system is unstable, for r = 10 years the 
system is stable, but passes through a phase of 
super-Eddington accretion, and for r = 0.1 year 
the system is stable and sub-Eddington. 

For weak tidal coupling (r = 10 15 years), shown 
in black, mass transfer causes the rotation rates of 
the donor and accretor to rapidly decrease and in¬ 
crease, respectively. For the case of stronger tidal 
coupling (r = 10 years and t = 0.1 year; shown in 
green and red, respectively), tidal forces exist to 
redistribute angular momentum, working to keep 
the spins of the donor and accretor synchronous 
with the orbit (/) = 1). 

The rapid increase in the semi-major axis for 



M a (M 0 ) 


Fig. 5.— As Figure 2, but using RL,asynch f° r the Roche 
lobe calculation. 


r = 10 15 years is a direct result of the conser¬ 
vation of angular momentum: the angular mo¬ 
mentum lost during the spin-down and mass loss 
of the donor is greater than the angular momen¬ 
tum gained by the accretor. The net decrease 
in spin angular momentum corresponds to an in¬ 
crease in the orbital angular momentum, increas¬ 
ing the semi-major axis. The mass transfer rate 
increases rapidly as the mass loss from the donor 
causes the radius to increase, even as the Roche 
lobe grows due to the increasing semi-major axis 
(see equations 34 and [35]). 


For the case of strong tidal coupling (r = 
10,0.1 year) we see a significant decrease in the 
mass transfer rate, with a correspondingly slower 
increase in the semi-major axis. The rotation rates 
of the component white dwarfs remain closer to 
synchronous, and we do not see the rapid runaway 
that is evident in the case of weak tidal coupling. 

Finally, we note that oscillations in the orbital 
parameters are seen in the r = 10 year case, which 
have been observed before by Gokhale2007. These 
oscillations are sensitive to small changes in the 
initial orbital parameters. We discuss this phe¬ 
nomenon further in section 4.5. 


4.4. Evolution of systems including the ef¬ 
fect of asynchronism on the Roche- 
lobe size 

Here, we show the results of our numerical in¬ 
tegrations using RL,asynch to calculate the Roche- 
lobe size. As discussed in section 2.6, the Eggleton 
approximation is dependent upon only the mass 
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Fig. 6.— As Figure 3, but with r = 10 years and using 

Rl ,asynch • 


ratio of the binary, whereas R.L,asynch accounts for 
deviations from synchronism. 

Figure 5 shows the end results of the systems 
with initial synchronization time-scales of r = 
10 15 years and evolved using RL,asynch- The color 
code remains the same as in Figures 2 and 3 and 
the dotted and dashed lines from Marsh2004 and 
Gokhale2007, respectively, for the 10 15 year time- 
scale are included for reference. As can be seen, 
this plot is nearly identical to Figure 2, which was 
calculated using the Eggleton approximation. 


Figure 6 shows the end results of the systems 
which begin with synchronization time-scales of 
r = 10 years. By comparing with Figure 3, 
we observe that, unlike the 10 15 year time-scale, 
Rl, asynch has a significant effect upon the end- 
states of the systems for the 10 year time-scale. In 
particular, the number of unstable black systems 
is reduced substantially. The cause of this differ¬ 
ence is explained in details below (section 4.5). 

In Figure 6, we also observe a slight increase 
in the number of stable systems across the entire 
range of accretor masses compared to Figure 3. 
Additionally, there is a small group of stable red 
systems beyond the widely-used q > 2/3 boundary 
for instability (for q > 2/3, the donor radius will 
expand faster than the Roche lobe which suggests, 
simplistically, that such systems will become un¬ 
stable). However, it is possible for systems with an 
initial mass ratio higher than 2/3 to survive, pro¬ 
vided that the mass ratio is close to the boundary 
(see, for example, D’Souza et al. ( |2006 >). Our 
results, as shown in Figure 6, confirm this obser- 



0 1000 2000 0 1000 2000 


time (years) 

Fig. 7. — Evolution of orbital parameters for r = 10 15 
years with initial masses of Md = 0.3 Mg and Mg = 
0.8 Mq. Solutions obtained using the Rl ,Egg calculation 
of the Roche lobe on left and the RL,asynch on the right. 


vation. 

/./. 1. Comparison of the evolution of systems 
using Rl ,Egg and RL,asynch solutions for 
t = 10 15 years 

Figure 7 compares the evolution of a system 
with an initial donor mass of 0.3 Mq and an ini¬ 
tial accretor mass of 0.8 Mq for r = 10 15 years. 
On left we show the evolution of various orbital 
parameters obtained using the Eggleton approxi¬ 
mation for the Roche lobe calculation. The right 
side mirrors the left but shows the evolution using 
Rl, asynch instead of Rl, Egg- As seen on the plots, 
the evolution of the orbital parameters happens on 
a much shorter time-scale for the Rl,E gg solution 
than for the RL,asynch solution. This is due to the 
fact that, as mass is transferred from the donor 
to the accretor, the rotation rate of the donor de¬ 
creases relative to the orbit and the rotation rate 
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Fig. 8.— Evolution of orbital parameters for r = 10 
years for initial masses, Mp = 0.38 Mq and Ma =0.8 Mq. 
First 1000 years of evolution with solutions obtained using 
^L,Egg f° r the calculation of the Roche lobe on left and 
using RL,asynch on right. The red lines show the super- 
Eddington accretion rate. 


of the accretor increases. As the donor spin ( f D ) 
decreases, A (which is proportional to equa¬ 
tion 39) also decreases. When using R,L,asynch, the 
size of the Roche lobe is inversely proportional to 
A. So, as A decreases, the Roche lobe will in¬ 
crease, which ultimately reduces the mass trans¬ 
fer rate. However, in the Rl ,Egg calculation, the 
Roche lobe is insensitive to changes in fn, so the 
mass transfer rate is higher relative to that cal¬ 
culated in the solution using RL,asynch- Since the 
mass transfer rate is higher in the Eggleton case, 
the orbital parameters change more quickly than 
in the RL,asynch case, as seen in Figure 7. 


4-4-2. Comparison of the evolution of systems 
using Rl,E gg and RL,asynch solutions for 
t = 10 years 

Figure 8 compares the evolution of a system 
with an initial donor mass of 0.38 Mq and an ini¬ 
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Fig. 9.— Full evolution in time of orbital parameter 
for r = 10 years with M g = 0.38 Mg,and Mg = 0.8 Mq. 
Compare to Figure 8 which shows only the first 1000 years 
of evolution for this particular system. 


tial accretor mass of 0.8 Mq with r = 10 years for 
the first 1000 years of evolution. Figure 9 shows 
the full evolution of this system, which evolved out 
to approximately 6 x 10 4 years before reaching a 
phase of disc accretion. As in Figure 7, the left 
panel shows the evolution using Rl,E gg and the 
right panel shows the evolution using RL,asynch- 
As can be seen in Figures 3 and 5, this system is 
stable under the RL,asynch calculation and super- 
Eddington under the Rl,E gg calculation. Unlike 
the 10 15 year time-scale, for the 10 year time- 
scale, the orbital parameters follow approximately 
the same path in both the Rl,E gg and RL,asynch 
cases, with the exception being the oscillations 
seen when using RL,E gg - We discuss these oscilla¬ 
tions in more detail below. 

4.5. Oscillations of orbital parameters 

As seen in Figure 8, using RL,asynch to calculate 
the Roche lobe damps the oscillations observed 
when using the Rl, Egg- In both cases, the mass 
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Fig. 10.— Evolution of the orbital period, P or b, com¬ 
pared to its initial orbital period, P or b,l , for a system with 
Mi) = 0.38 Mg and Mg = 0.8 Mg for various synchroniza¬ 
tion time-scales at contact. Red, blue, and magenta lines 
represent the orbital periods for synchronization timescales 
of 20, 30, and 40 years, respectively. The left-hand panel 
shows the evolution using Rz,,Egg and the right-hand panel 
shows the evolution using RL,asynch- 


transfer rate increases initially, which causes the 
accretor to spin up and the donor to spin down. 
In the Rl, asynch case, this causes the mass trans¬ 
fer rate to increase at a slower rate relative to the 
Rl, Egg case. [^] As discussed above, when using 
Rl Eggj the mass transfer rate will continue to in¬ 
crease until tides have transferred sufficient angu¬ 
lar momentum from the spin of the components to 
the orbit to widen the orbit and reduce the mass 
transfer rate. Once this occurs, the mass transfer 
rate is reduced and tides continue to redistribute 
angular momentum between the component spins 
and the orbit (in this case the donor is spun up 
while the accretor is spun down). But again, be¬ 
cause changes in the donor spin are not accounted 
for when using RL,Egg, the mass transfer rate will 
again “over shoot” and the oscillations continue. 
The dependence of the Roche lobe radius upon the 
donor spin when using RL,asynch slows the changes 
in the mass transfer rate and damps the oscilla¬ 
tions. 

Although Figure 8 shows that the short-term 
evolution of this particular system differs substan¬ 
tially between the Rl, Egg and RL,asynch Roche 
lobe calculations, Figure 9 shows that the long¬ 
term evolution of this system is very similar. Al¬ 
though the oscillations do not impact the long¬ 
term evolution of many of the systems, the oscil¬ 
lations do impact whether a system is categorized 


5 As described previously, as fo decreases, A decreases, 
which causes the Roche lobe to increase and reduces the 
mass transfer rate. 


PL, 


PL, 


PL 


CL 



Fig. 11.— Evolution of the orbital period for two dif¬ 
ferent mass ratios. The upper graphs show the evolution 
for an initial mass ratio of q = 0.425 (Mo = 0.34 Mg and 
Mjl i = 0.8 Mg) and the lower graphs show an initial mass 
ratio of q=0.475 (Mo = 0.38 Mq and Ma = 0.8 Mg). 
The left-hand panel shows evolution using Ro,Egg and the 
righ-hand panel shows evolution using RL,asynch■ As in 
Figure 10, the red plots show evolution for an initial syn¬ 
chronization time-scale of r = 10 yrs and the magenta plots 
show the evolution for r = 40 yrs. 


as sub-Eddington or super-Eddington or whether 
it reaches unstable levels of mass transfer. In the 
case of the system illustrated in Figures 8 and 9, 
the oscillations cause the system briefly to become 
super-Eddington, while the system would remain 
sub-Eddington in the absence of oscillations. 

In general, as the mass transfer rates spike 
initially under the Rl ,Egg calculation before ul¬ 
timately settling to the RL,asynch value, it may 
reach super-Eddington or unstable levels, while, in 
the RL,asynch case, the mass transfer rate may re¬ 
main stable instead of super-Eddington, or super- 
Eddington instead of unstable. This explains the 
lack of unstable black systems in Figure 5 com¬ 
pared to Figure 3. For the systems shown in Fig¬ 
ure 3, the oscillations in the orbital parameters 
arising from the use of RL,Egg may cause the mass 
transfer rate to reach artificially high values before 
the oscillations relax, leading to the large popu¬ 
lation of black systems at the top of the figure. 
Meanwhile, since using RL, aS ynch damps the oscil¬ 
lations, the mass transfer rate is able to stay below 
the limit for stability of 0.01 Mq/ yr. 

Both Marsh2004 and Gokhale2007 noted the 
presence of oscillations in the orbital parame¬ 
ters near the boundary between stable and un- 


16 








stable systems. These analyses observed that, for 
such systems, small changes in the synchroniza¬ 
tion time-scales and in the mass ratio will cause 
changes in the frequency and amplitude of the 
oscillations. Figures 10 and 11 show the evolu¬ 
tion of the orbital period for various synchroniza¬ 
tion time-scales and mass ratios, respectively, us¬ 
ing both Rl,E gg and RL,asynch to calculate the 
size of the Roche lobe. These figures confirm the 
observations of Marsh2004 and Gokhale2007 that 
the amplitude and frequency of the oscillations is 
dependent upon both the mass ratio and the syn¬ 
chronization time-scale. 

The left-hand panel of Figure 10 illustrates that 
the amplitude of the oscillations increases with in¬ 
creasing initial synchronization time-scale, while 
the frequency of the oscillations decreases. These 
oscillations can be explained as follows: Once mass 
transfer begins, if the system is unstable or nearly 
unstable, the time-scale for changes in the mass 
transfer rate becomes greater than the time-scale 
at which the semi-major axis increases due to 
tides. As a result, Md increases rapidly and the 
accretor spin increases significantly. The result is 
a large asynchronicity between the accretor, the 
orbit, and the donor. As the orbital separation 
increases due to the effect of tides, Mo decreases, 
until the timescale for changes in the mass transfer 
rate becomes less than the time-scale at which the 
semi-major axis increases due to tides. At this 
point, tides are able to redistribute angular mo¬ 
mentum from the spin of the accretor back into the 
orbit. If enough angular momentum is transferred 
back into the orbit, the donor becomes detached. 
Once tides effectively re-synchronize the spins of 
the stars with the orbit, gravitational wave losses 
act to bring the binary back into contact, and the 
entire cycle repeats. 

The right-hand panel of Figure 10 demonstrates 
that, for this particular mass ratio, the oscilla¬ 
tions are not present when using RL,asynch to cal¬ 
culate the size of the Roche lobe. Furthermore, 
the right-hand panel illustrates that, in the ab¬ 
sence of oscillations, the orbital period, and there¬ 
fore, the semi-major axis increases more rapidly 
for lower r values of the synchronization timescale. 
For lower r values, the tidal coupling is stronger, 
which means spin angular momentum is able to 
be re-distributed more efficiently, which allows the 
semi-major axis to remain larger compared to the 


semi-major axis of systems evolved using higher 
values of r. 

Figure 11 demonstrates the effect of mass ra¬ 
tio upon the amplitude of the oscillations. As 
Gokhale2007 observes, we see that for a fixed r 
value, a lower mass ratio (top two graphs of Fig¬ 
ure 11) yields oscillations of lower amplitude com¬ 
pared to a higher mass ratio (lower two graphs of 
Figure 11). In Figure 3, we see that switching from 
a mass ratio of 0.475 to 0.425 moves us away from 
the boundary which separates super-Eddington 
systems from stable systems. This demonstrates 
that, in general, as we move away from the stabil¬ 
ity boundary, the oscillations are reduced. 

The two right-hand panels of Figure 11 show 
the evolution using RL,asynch for the roche lobe. 
This again demonstrates that oscillations are not 
present when the Roche lobe depends on the asyn¬ 
chronicity between the donor and orbit. 


5. Conclusions 


We have studied the long-term evolution of 
DWD binary systems undergoing direct-impact 
mass transfer, including the effects due to mass 
transfer, gravitational radiation, and tidal forces 
arising from asynchronicity between the donor, ac¬ 
cretor, and the orbit. We implemented the ballis¬ 


tic mass-transfer treatment developed in (Sepin- 


sky et al. 2010) to calculate the changes to or¬ 


bital parameters during direct-impact mass trans¬ 
fer. By implementing this method, we found that 
the number of stable DWD systems increased for 
both weak and strong tidal coupling compared 
with the results of Marsh2004 and Gokhale2007, 
as shown in Figures 2 and 3. 


For the first time, we also account for the mod¬ 
ification of the Roche-lobe size due to the asyn¬ 
chronicity of the donor. As a result, we find that 
the number of stable systems increases, particu¬ 
larly for the case of strong tidal coupling, as shown 
in Figures 5 and 6. When not accounting for the 
asynchronicity effects on the Roche-lobe size, we 
reproduce the oscillations in the orbital parame¬ 
ters first noted by Goklrale2007. We found that, 
when the size of the Roche lobe was permitted to 
vary with donor asynchronicity, these oscillations 
were dampened, as shown in Figure 8. We con¬ 
clude that the oscillations created when using the 
Eggleton approximation for the Roche lobe calcu- 
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lation create artificially high mass transfer rates 
which leads to an artificially high number of super- 
Eddington and unstable systems. By eliminat¬ 
ing the oscillations, our treatment yields a higher 
number of sub-Eddington and stable systems. 


We expect systems which are stable through¬ 
out our calculations (red systems in Figures 2, 3, 
5, and 6) to be AM CVn progenitors. As a re¬ 
sult of the increase in stable systems shown here 
compared to previous analyses, we conclude that 
DWD evolution may be a more likely avenue for 
the creation of AM CVn than previously expected. 


In future analyses, we intend to investigate 
the case of initially eccentric binaries with asyn¬ 
chronous component stars, as we expect there may 
be a significant population of DWDs that have 
not had time to circularize and synchronize by the 


time mass transfer occurs (Willems et al. 2007). 


Additionally, we intend to implement a treatment 
of disc accretion so that we can track the evolution 
of the DWD systems through all types of accretion 
processes. 
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